///////////////////////////////////////////////////////////////////////////////
/******************************************************************************

qdd: A Quantum Computer C++ Library constructed using BDD's

Copyright (C) 1999-2007 David Greve <TheBeaNerd@gmail.com>

This library is free software; you can redistribute it and/or modify
it under the terms of the GNU Library General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.

This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Library General Public License for more details.

You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA.

******************************************************************************/
///////////////////////////////////////////////////////////////////////////////

#include "qdd.h"

QuantumPhase::QuantumPhase

///////////////////////////////////////////////////////////////////////////////
/******************************************************************************

  class::QuantumException

******************************************************************************/
///////////////////////////////////////////////////////////////////////////////

QuantumException::QuantumException (const char *string) {

  msg = string;

}

/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  class::QuantumCondition

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

QuantumCondition::QuantumCondition () {

  QDD_ECHO ( cout << "QuantumCondition()\n" ) ;

  ConditionVector = bdd_false();

}

QuantumCondition & QuantumCondition::operator & (const QuantumBit & arg) {

  QDD_ECHO ( cout << "cond = cond & qbit\n" ) ;

  ConditionVector = bdd_ithvar(arg.Index);

  return(* this) ;

}

QuantumCondition & QuantumCondition::operator & (QuantumCondition & arg) {

  QDD_ECHO ( cout << "cond = cond & cond\n" ) ;

  ConditionVector |= arg.ConditionVector

  return(* this);

}


/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  class::QuantumBit

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

double QuantumBit::Denominator = 1.0;
bdd  QuantumBit::State = bddfalse;
int *QuantumBit::IndexHeap;
int  QuantumBit::NextFreeIndex = -1;

int  QuantumBit::INITIALIZED = 0;
int  QuantumBit::MIN_VAR = -1;
int  QuantumBit::MAX_VAR = -1;

// DAG : Something here ain't right .. do we need to declare "friend"
// functions outside of the class description?

extern Integer Ipow(long, long);

void qdd_init(int QBits, int NBits, int megs) {

  qdd_assert(! QuantumBit::INITIALIZED, "qdd_init: Redundant QDD initialization");

  int QDD_BDD_NODENUM = (int) (megs * (2000000.0/64.0)) ;

  int QDD_BDD_CACHESIZE = (int) (QDD_BDD_NODENUM / 10.0) ;

  // So .. we need a var for complex numbers.
  // We need some number of vars for the mantissa
  // We need a var for mantissa overflow
  // We need one var per quantum bit (which includes any temporary bits).

  QuantumBit::NBits = NBits;
  QuantumBit::PBits   = ilog2(2*NBits + 1);
  QuantumBit::SignBit = ipow2(QuantumBit::PBits) - 1;
  QuantumBit::ZVar  = PBits;

  // Total number of auxillary bits needed.
  int AuxBitCount = PBits + 1;

  QuantumBit::MIN_VAR = AuxBitCount;
  QuantumBit::MAX_VAR = AuxBitCount + QBits;

  if (QDD_BDD_NODENUM < (Ipow(2,MAX_VAR) * (int) (2000000.0/(2 << 16)))) {

    cerr << "** Warning : RAM to QBIT ratio dangerously small\n" ;

  }

  bdd_init(QDD_BDD_NODENUM,QDD_BDD_CACHESIZE);
  bdd_setvarnum(QuantumBit::MAX_VAR);
  bdd_autoreorder(BDD_REORDER_WIN2);
  bdd_reorder_verbose(1);

  QuantumBit::IndexHeap = new int[QuantumBit::MAX_VAR];

  int i;
  for (i=0;i<QuantumBit::MAX_VAR;i++) {
    QuantumBit::IndexHeap[i] = i;
  }
  QuantumBit::NextFreeIndex = QuantumBit::MIN_VAR;
  
  srand(7);
  
  QuantumBit::State = Load(QuantumBit::PBits,
			   QuantumBit::SignBit,
			   QuantumBit::NBits,
			   QuantumBit::MAX_VAR,
			   0,
			   (int)exp2(QuantumBit::NBits-1),
			   QuantumBit::State);

  QuantumBit::Denominator = 1.0/exp2(QuantumBit::NBits-1);

  QuantumBit::INITIALIZED = 1;

}

void QuantumBit::QuantumBitInit(int value) {

  QDD_ECHO ( cout << "QuantumBitInit(" << value << ")\n" ) ;

  qdd_assert(((value == 0) | (value == 1)), "QuantumBitInit: Initial QBit value must be 0 or 1");

  qdd_assert(INITIALIZED, "QuantumBitInit: Uninitialized Quantum State");

  qdd_assert((NextFreeIndex < MAX_QBITS),"QuantumBitInit: Exhausted Quantum Bits");

  int index = IndexHeap[NextFreeIndex++];

  Index = index;

  Set(value);

}

QuantumBit::QuantumBit() {

  QDD_ECHO ( cout << "QuantumBit()\n" ) ;

  QuantumBitInit(0);

}

QuantumBit::QuantumBit(int value) {

  QDD_ECHO ( cout << "QuantumBit(" << value << ")\n" ) ;

  QuantumBitInit(value);

}

QuantumBit::~QuantumBit() {

  QDD_ECHO ( cout << "qbit[" << Index << "].~QuantumBit()\n" ) ;

  // Perform Quantum Garbage Collection

  if (! this.Pure()) {
    cerr << "** Warning : ~qbit() : Collecting Mix()'d Bit\n" ;
    Sample();
  }

  putValue(bddfalse);

  qdd_assert(NextFreeIndex > MIN_USER_QBIT,
	 "~qbit(): Excessive DeAllocation (?)");

  qdd_assert(Index >= MIN_USER_QBIT,
	 "~qbit() : Corrupt Quantum Bit");

  IndexHeap[--NextFreeIndex] = Index;

  Index = -1;

}

// DAG - eliminated
// void QuantumBit::operator ~ ()

void QuantumBit::Not () {
  bdd x  = bdd_ithvar(Index);
  bdd Tx = bdd_restrict(State,x);
  bdd Fx = bdd_restrict(State,!x);
  State = bdd_ite(x,Fx,Tx);
}

void QuantumBit::operator ^= (QuantumCondition &arg) {

  QDD_ECHO ( cout << "qbit[" << Index << "] ^= (cond)\n" ) ;

  bdd x = bdd_ithvar(Index);
  bdd C  = bdd_restrict(State,arg.ConditionVector);
  bdd Tx = bdd_restrict(State,x);
  bdd Fx = bdd_restrict(State,!x);
  State = bdd_ite(x ^ C,Tx,Fx);

}

void QuantumBit::operator ^= (const QuantumBit &arg) {

  QDD_ECHO ( cout << "qbit[" << Index << "] ^= (qbit[" << arg.Index << "])\n" ) ;

  bdd x = bdd_ithvar(Index);
  bdd ConditionVector = bdd_ithvar(arg.Index);
  bdd C  = bdd_restrict(State,arg.ConditionVector);
  bdd Tx = bdd_restrict(State,x);
  bdd Fx = bdd_restrict(State,!x);
  State = bdd_ite(x ^ C,Tx,Fx);

}

// DAG - eliminated
// QuantumCondition & QuantumBit::operator ! ()
// QuantumCondition & QuantumBit::operator & (const QuantumBit & arg)
// QuantumCondition & QuantumBit::operator & (QuantumCondition & arg)

void QuantumBit::Mix() {

  QDD_ECHO ( cout << "qbit[" << Index << "].Mix()\n" ) ;

  bdd x = bdd_ithvar(Index);
  bdd Tx = bdd_restrict(State,x);
  bdd Fx = bdd_restrict(State,!x);

  // Let's make phase a class unto itself.  This class can contain the
  // Denominator as well as the boundaries on the phase such as pbits,
  // Nbits, etc.

  // The Index should be stored as a bdd.  We can do checks on "unit"
  // bdds in appropriate locations.

  Sx = bdd_add(Tx,Fx);
  Dx = bdd_sub(Tx,Fx);

  Denominator *= sqrt(2.0);

  State = bdd_ite(x,Dx,Sx);
  

}

bool QuantumBit::Pure() {

  bdd x = bdd_ithvar(Index);

  return((bdd_restrict(State, x) == bddfalse) ||
	 (bdd_restrict(State,!x) == bddfalse))

}

int QuantumBit::Sample() {

  bdd x = bdd_ithvar(Index);

  double p = PChance(PBits,NBits,x,State);

  double r = (1.0 * rand()) / RAND_MAX;

  int value = 0;

  if (< r p) {
    State = bdd_restrict(State,x);
    Denominator *= 1.0-p;
    value = 1;
  } else {
    State = bdd_restrict(State,!x);
    Denominator *= p;
    value = 0;
  }

  return(value);

}

void QuantumBit::Set(int value) {

  bdd Tx = bdd_restrict(State,bdd_ithvar(Index));
  bdd Fx = bdd_restrict(State,bdd_nthvar(Index));

  QDD_ECHO ( cout << "qbit[" << Index << "].Set(" << value << ")\n" ) ;

  qdd_assert((value == 1) || (value == 0), "qbit.Set(int) takes a boolean value");

  if (Tx == bdd_false()) {
    State = value ? bdd_ite(x,Fx,Tx) : State;
  } else if (Fx == bdd_false()) {
    State = value ? State : bdd_ite(x,Fx,Tx);
  } else {
    qdd_assert(0 == 1,"QuantumBit::Set(): Attempting to set an Impure state");
  }

}

/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  class::QuantumRegister

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

QuantumRegister::QuantumRegister(int size) {

  QDD_ECHO ( cout << "QuantumRegister( " << size << " )\n" ) ;

  qdd_assert( size > 0, "qreg(): Register Size must be greater than zero");

  QuBit = new QuantumBit[size];

  Size = size;

}

QuantumRegister::~QuantumRegister() {

  QDD_ECHO ( cout << "~QuantumRegister()\n" ) ;

  delete [] this->QuBit;

}

QuantumBit & QuantumRegister::operator [] (int index) {

  qdd_assert(((index >= 0) && (index < Size)),"qreg[] : Index out of range");

  return(QuBit[index]) ;

}

Integer QuantumRegister::Sample() {

  QDD_ECHO ( cout << "qreg.Sample()\n" ) ;

  Integer value = 0;

  int i;
  for (i=(Size-1);i>=0;i--) {
    value <<= 1;
    value += QuBit[i].Sample();
  }

  return(value);

}

void QuantumRegister::Set(const Integer & value) {

  QDD_ECHO ( cout << "qreg.Set(" << value << ")\n" ) ;

  int i;
  for (i=0;i<Size;i++) {
    QuBit[i].Set(BIT(value,i)) ;
  }

}

void QuantumRegister::operator ^= (const QuantumRegister & arg) {

  QDD_ECHO ( cout << "qreg ^= qreg\n" ) ;

  int i;
  for (i=0;i<Size;i++) {
    QuBit[i] ^= arg.QuBit[i];
  }

}

void QuantumRegister::Mix() {

  QDD_ECHO ( cout << "qreg.Mix()\n" ) ;

  int i;
  for (i=0;i<Size;i++) {
    QuBit[i].Mix();
  }

}

void QuantumRegister::Pure() {

  QDD_ECHO ( cout << "qreg.Pure()\n" ) ;

  int i;
  for (i=0;i<Size;i++) {
    QuBit[i].Pure();
  }

}

Integer QuantumRegister::Period() {

  QDD_ECHO ( cout << "qreg.Period()\n" ) ;

  bdd s0 = bddtrue;
  Integer v0 = 0;
  int i;
  for (i=(Size-1);i>=0;i--) {
    v0 <<= 1;
    if (bdd_simplify(QuBit[i].Value(),s0) == bddtrue) {
      v0 += 1;
      s0 &= QuBit[i].Value(); // I don't think we need this.
    }
    else {
      s0 &= (! QuBit[i].Value());
    }
  }

  bdd s1 = (! s0);
  Integer v1 = 0;
  for (i=(Size-1);i>=0;i--) {
    v1 <<= 1;
    if (bdd_simplify(QuBit[i].Value(),s1) == bddtrue) {
      v1 += 1;
      s1 &= QuBit[i].Value();
    }
    else {
      s1 &= (! QuBit[i].Value());
    }
  }  

  (* this).Sample();

  Integer period = v1 - v0;

  if (s1 == bddfalse) {
    period = Ipow(2,Size);
  }

  return(period);

}
